{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brief Introduction to Hyperspectral Unmixing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of hyperspectral unmixing is to decompose an image cube into the pure spectral signatures found in a scene (called endmembers) and the amount (or proportion) of each endmember found in each pixel. This is sub-pixel analysis since we are determining sub-pixel amounts of each material in each pixel.\n",
    "\n",
    "When performing hyperspectral unmixing, we first must assume a particular mixing model.  \n",
    "\n",
    "The most common mixing model used in practice is the *Linear Mixing Model* (also known as the *Convex Geometry Model*).  Although it is the most commonly used, it often does not hold in practice.  \n",
    " \n",
    "<img src=\"Picture3.png\" alt=\"Hyperspectral Mixing Models\" style=\"width: 700px;\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a number of non-linear mixing models to account for canopies and multi-level mixing and intimate mixing in imagery. These models include: \n",
    "<ul> \n",
    "<li> *Hapke, Kulbelka-Munk and Shkuratov Models*: Physics-based mixing models relying on radiative transfer theory.  Computationally complex and requires significant knowledge of scene parameters to perform accurately. \n",
    "<ul>\n",
    "<li> R. Close, P. Gader, J. Wilson, A. Zare, \"Using physics-based macroscopic and microscopic mixture models for hyperspectral pixel unmixing\", Proc. SPIE 8390, Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVIII, 83901L (24 May 2012); doi: 10.1117/12.919583; <url> http://dx.doi.org/10.1117/12.919583</url>\n",
    "<li> B. Hapke, “Bidirection reflectance spectroscopy. I. theory,” J. Geo- phys. Res., vol. 86, pp. 3039–3054, 1981.\n",
    "<li> P. Kulbelka and F. Munk, “Reflection characteristics of paints,”\n",
    "Zeitschrift fur Technische Physik, vol. 12, pp. 593–601, 1931.\n",
    "<li> Y. Shkuratov, L. Starukhina, H. Hoffmann, and G. Arnold, “A model of spectral albedo of particulate surfaces: Implications for optical properties of the Moon,” Icarus, vol. 137, p. 235246, 1999.\n",
    "</ul>\n",
    "<li> *Piece-wise Convex Mixing*: Represent scene with discrete sets of linear mixtures.  Accounts for disparate regions in scene (e.g., an image covering urban and rural regions will likely have two distinct sets of endmembers associated with each region).  \n",
    "<ul>\n",
    "<li> A. Zare, P. Gader, O. Bchir and H. Frigui, \"Piecewise Convex Multiple-Model Endmember Detection and Spectral Unmixing,\" in IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 5, pp. 2853-2862, May 2013. <url>http://ieeexplore.ieee.org/abstract/document/6352892/</url>\n",
    "<li> A. Zare, O. Bchir, H. Frigui and P. Gader, \"Spatially-smooth piece-wise convex endmember detection,\" 2010 2nd Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing, Reykjavik, 2010, pp. 1-4. <url>http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=5594897&isnumber=5594823</url> \n",
    "</ul>\n",
    "<li> *Non-physics/Manifold Based*: Represent non-linearities in data with non-linear models commonly used in statistical machine learning literature such as kernel approaches, non-linear manifold learning and others. \n",
    "<ul>\n",
    "<li> K. J. Guilfoyle M. L. Althouse C.-I. Chang \"A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks\" IEEE Trans. Geosci. Remote Sensing, vol. 39 no. 8 pp. 2314-2318 Aug. 2001. <url>http://ieeexplore.ieee.org/document/957296/</url>\n",
    "<li> A. Halimi, Y. Altmann, N. Dobigeon and J. Y. Tourneret, \"Nonlinear Unmixing of Hyperspectral Images Using a Generalized Bilinear Model,\" in IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 11, pp. 4153-4162, Nov. 2011. <url> http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=5702384&isnumber=6059695</url>\n",
    "<li> Y. Altmann N. Dobigeon S. McLaughlin J.-Y. Tourneret \"Nonlinear unmixing of hyperspectral images using radial basis functions and orthogonal least squares\" Proc. IEEE Int. Conf. Geoscience and Remote Sensing (IGARSS) pp. 1151-1154 July 2011. <url>http://ieeexplore.ieee.org/document/6049401/</url>\n",
    "<li> P. Gader D. Dranishnikov A. Zare J. Chanussot \"A sparsity promoting bilinear unmixing model\" Proc. IEEE GRSS Workshop Hyperspectral Image Signal Processing: Evolution Remote Sensing (WHISPERS), June 2012. <url>http://ieeexplore.ieee.org/document/6874255/</url>\n",
    "<li> and many others..\n",
    "</ul>\n",
    "<li> *Overview of non-linear mixing*: \n",
    "<ul>\n",
    "<li>R. Heylen, M. Parente and P. Gader, \"A Review of Nonlinear Hyperspectral Unmixing Methods,\" in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 7, no. 6, pp. 1844-1868, June 2014. <url>http://ieeexplore.ieee.org/abstract/document/6816071/</url>\n",
    "<li> N. Dobigeon, J. Y. Tourneret, C. Richard, J. C. M. Bermudez, S. McLaughlin and A. O. Hero, \"Nonlinear Unmixing of Hyperspectral Images: Models and Algorithms,\" in IEEE Signal Processing Magazine, vol. 31, no. 1, pp. 82-94, Jan. 2014. <url>http://ieeexplore.ieee.org/abstract/document/6678284/</url>\n",
    "</ul>\n",
    "</ul>\n",
    "\n",
    "In addition to non-linear mixing, the linear mixing model may not hold when considering spectral variability.   Spectral variability can be caused by environmental conditions (e.g., variations in illumination), atmospheric conditions (e.g., water in atmosphere), and inherent variability within a material.  Inherent variability depends on the scale of the endmember under consideration. For example, if a particular plant species is associated to one endmember, variation in this endmember may occur due to the upper and under-side of leaves of that species having different spectral signatures).  Spectral unmixing methods that account for spectral variability can be organized into two categories: set-based approaches and distribution-based approaches.  Set-based approaches represent an endmember using a discrete set of endmember spectra. Distribution-based approaches use a full probability distribution to represent an endmember and its associated variability.  Often, set-based approaches under-represent the variability whereas distribution-based approaches may over-represent the variability.   Examples of unmixing methods that account for spectral variability include: \n",
    "<ul>\n",
    "<li> *MESMA*: A set-based approach, Multiple Endmember Spectral Mixture Analysis, \n",
    "<li> *AAM*: A set-based approach, Alternating Angle Minimization:  R. Heylen, A. Zare, P. Gader and P. Scheunders, \"Hyperspectral unmixing with endmember variability via alternating angle minimization,\" IEEE Tran. Geosci. Remote Sens., vol. 54, no. 8, pp. 4983-4993, Aug. 2016. Paper: <url>http://ieeexplore.ieee.org/document/7464927/</url> Code: <url>https://sites.google.com/site/robheylenresearch/code/AAM.zip?attredirects=0&d=1</url>\n",
    "<li> *Normal Compositional Model*: A distribution-based approach where each endmember is represented using a Gaussian distribution.  There are a number of algorithms based on the NCM including: \n",
    "<ul>\n",
    "<li> D. Stein, \"Application of the normal compositional model to the analysis of hyperspectral imagery,\" IEEE Workshop on Advances in Techniques for Analysis of Remotely Sensed Data, 2003, 2003, pp. 44-51. <url> http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1295171&isnumber=28800</url>\n",
    "<li> O. Eches, N. Dobigeon, C. Mailhes and J. Y. Tourneret, \"Bayesian Estimation of Linear Mixtures Using the Normal Compositional Model. Application to Hyperspectral Imagery,\" in IEEE Transactions on Image Processing, vol. 19, no. 6, pp. 1403-1413, June 2010. <url>http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=5427031&isnumber=5464460</url>\n",
    "<li> A. Zare, P. Gader and G. Casella, \"Sampling Piecewise Convex Unmixing and Endmember Extraction,\" in IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 3, pp. 1655-1665, March 2013. <url>http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6297456&isnumber=6469260</url>\n",
    "</ul>\n",
    "<li> *Beta Compositional Model*: A distribution-based approach where each endmember (and each band/wavelength) is represented using a Beta distribution to enforce endmember reflectance values remain between 0 and 1.   Paper: X. Du, A. Zare, P. Gader and D. Dranishnikov, \"Spatial and Spectral Unmixing Using the Beta Compositional Model,\" in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 7, no. 6, pp. 1994-2003, June 2014. <url>http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6851850&isnumber=6870503</url>\n",
    "<li> *Overview papers on unmixing given spectral variability*: \n",
    "<ul>\n",
    "<li>A. Zare and K. C. Ho, \"Endmember Variability in Hyperspectral Analysis: Addressing Spectral Variability During Spectral Unmixing,\" in IEEE Signal Processing Magazine, vol. 31, no. 1, pp. 95-104, Jan. 2014. <url>http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6678271&isnumber=6678221</url>\n",
    "<li> Somers, Ben, et al. \"Endmember variability in spectral mixture analysis: A review.\" Remote Sensing of Environment 115.7 (2011): 1603-1616. <url>https://www.sciencedirect.com/science/article/pii/S0034425711000800</url>\n",
    "</ul>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear mixing model assumes each pixel is a convex combination of pure *endmember* spectra.   Endmembers are the spectral signatures of the pure, constituent materials in a scene.  The linear mixing model can be written as: \n",
    "\n",
    "$\\mathbf{x}_i = \\sum_{k=1}^M p_{ik}\\mathbf{e}_{k} + \\epsilon_i \\quad i= 1, \\ldots, N$\n",
    "\n",
    "where $N$ is the number of pixels in the image, $M$ is the number of endmembers, $\\epsilon_i$ is the residual error term, $p_{ik}$ is the *proportion* (also called *abundance*) of the $k$th endmember in the $i$th pixel, $\\mathbf{e}_k$ is the spectral signature of the $k$th endmember, and $\\mathbf{x}_i$ is the spectral signature of the $i$th pixel. \n",
    "\n",
    "In this model, the proportions are assumed to sum to one and be non-negative (as they refer to percentages of material found within a pixel): \n",
    "\n",
    "$p_{ik} \\ge 0 \\quad \\forall i,k$\n",
    "\n",
    "$\\sum_{k=1}^M p_{ik} = 1$\n",
    "\n",
    "The linear mixing model (also sometimes called the \"Convex Geometry Model\" can be visualized as shown in the image below.  Under this model, each pixel lies within the convex hull defined by the endmembers.  Also, the endmembers are called *endmembers* because they are found out at the ends of the data.  It has been shown that this model is effective at modeling mixtures due to inadequate spatial resolution by the hyperspectral imager (but not due to mixing on the ground or multiple reflections). \n",
    "\n",
    "<img src=\"Picture04.png\" alt=\"Linear Mixing Model\" style=\"width: 400px;\"/>\n",
    "\n",
    "Due to the linear mixing model, we often have the goal of \"unmixing\" a hyperspectral data cube.  The goal in unmixing is to, given the data $\\mathbf{X} = \\left\\{ \\mathbf{x}_i \\right\\}_{i=1}^N$, estimate the endmember spectral signatures and their proportions founds within each pixel in a hyperspectral data cube.  Note, this problem amounts to an ill-posed matrix factorization problem.  Thus, to solve it, we generally have to impose constraints on the endmebmers and proportions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# imports and setup\n",
    "import numpy as np\n",
    "import os.path\n",
    "import scipy.io\n",
    "from loadmat import loadmat\n",
    "\n",
    "import matplotlib as mpl\n",
    "default_dpi = mpl.rcParamsDefault['figure.dpi']\n",
    "mpl.rcParams['figure.dpi'] = default_dpi*2\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load gulfport campus image\n",
    "img_fname = 'muufl_gulfport_campus_w_lidar_1.mat'\n",
    "spectra_fname = 'tgt_img_spectra.mat'\n",
    "\n",
    "dataset = loadmat(img_fname)['hsi']\n",
    "\n",
    "hsi = dataset['Data'][:,:,4:-4] # trim noisy bands \n",
    "valid_mask = dataset['valid_mask'].astype(bool)\n",
    "n_r,n_c,n_b = hsi.shape\n",
    "wvl = dataset['info']['wavelength'][4:-4]\n",
    "rgb = dataset['RGB']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After loading the data, lets extract some endmembers using the Pixel Purity Index algorithm.  This algorithm assumes that pure spectra for each endmember can be found in the scene.  This assumption that does not hold for highly mixed data sets. \n",
    "\n",
    "Reference for PPI: J. W. Boardman, \"Automated spectral unmixing of AVIRIS data using convex geometry concepts\", Summaries 4th JPL Airborne Geoscience Workshop, Jet Propulsion Lab., vol. 1, pp. 11-14, 1993.\n",
    "\n",
    "Of course, there are MANY algorithms in the literature besides PPI that estimate endmember spectra.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# extract some endmembers using Pixel Purity Index algorithm\n",
    "#  using PySptools from https://pysptools.sourceforge.io\n",
    "# Exercise: Vary number of endmembers and number of skewers to see changes in endmember estimation\n",
    "import pysptools\n",
    "import pysptools.eea\n",
    "\n",
    "hsi_array = np.reshape(hsi,(n_r*n_c,n_b))\n",
    "valid_array = np.reshape(valid_mask,(n_r*n_c,))\n",
    "M = hsi_array[valid_array,:]\n",
    "q = 5 #Number of Endmembers\n",
    "numSkewers = 500 #PPI parameter of number of projections used to find extreme data points that may be endmembers\n",
    "E,inds = pysptools.eea.eea.PPI(M, q, numSkewers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# plot the endmembers we found\n",
    "plt.plot(wvl,E.T)\n",
    "plt.xlabel('wavelength (nm)')\n",
    "plt.ylabel('reflectance')\n",
    "plt.legend([str(i+1) for i in range(q)])\n",
    "plt.title(\"PPI Endmembers\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After estimating endmember spectra, we can estimate the abundances/proportions for each pixel in the image.  We will use the FCLS algorithm for this.  (Again, there are many algorithms in the literature that estimate proportions given endmembers.  FCLS is just one example.)\n",
    "\n",
    "Reference for FCLS: D. C. Heinz and Chein-I-Chang, \"Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery,\" in IEEE Transactions on Geoscience and Remote Sensing, vol. 39, no. 3, pp. 529-545, Mar 2001.\n",
    "<url> http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=911111&isnumber=19663 </url>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# find abundances given the endmembers\n",
    "import pysptools.abundance_maps\n",
    "\n",
    "maps = pysptools.abundance_maps.amaps.FCLS(M, E) #This runs slowly with large data sets/more endmembers. \n",
    "#maps = np.zeros((M.shape[0],E.shape[1])) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# re-ravel abundance maps\n",
    "map_imgs = []\n",
    "for i in range(q):\n",
    "    map_lin = np.zeros((n_r*n_c,))\n",
    "    map_lin[valid_array] = maps[:,i]\n",
    "    map_imgs.append(np.reshape(map_lin,(n_r,n_c)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# display abundance maps\n",
    "for i in range(q):\n",
    "    plt.figure()\n",
    "    plt.imshow(map_imgs[i],vmin=0,vmax=1)\n",
    "    plt.colorbar()\n",
    "    plt.title('FCLS Abundance Map %d'%(i+1,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can estimate endmembers, number of endmembers and abundances simultaneously using the SPICE algorithm.  SPICE is also applicable to highly mixed datasets as it does not assume endmember spectra can be found within the data set.  Of course, this is only one example of this type of algorithm in literature. \n",
    "\n",
    "Reference for SPICE: Zare, A.; Gader, P.; , \"Sparsity Promoting Iterated Constrained Endmember Detection in Hyperspectral Imagery,\"\" IEEE Geoscience and Remote Sensing Letters, vol.4, no.3, pp.446-450, July 2007.\n",
    "    <url>https://faculty.eng.ufl.edu/machine-learning/2007/07/zare2007sparsitypromoting/</url>\n",
    "    \n",
    "Matlab code for SPICE can be found here: <url>https://github.com/GatorSense/SPICE</url>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# run SPICE to find number of endmembers, endmembers, and abundances simultaneously\n",
    "#Exercise: Vary SPICE parameters (in params) to see effect on endmember and parameter estimation. \n",
    "from SPICE import *\n",
    "\n",
    "params = SPICEParameters()\n",
    "inputData = M.T.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# to save time, downsample inputData\n",
    "dsData = inputData[:,::20]\n",
    "dsData.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# run SPICE\n",
    "[eM,dsP] = SPICE(dsData,params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# unmix endmembers again with full data matrix (because we downsampled for sake of time)\n",
    "P = unmix2(inputData,eM)\n",
    "n_em = eM.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#plot endmembers\n",
    "plt.plot(wvl,eM)\n",
    "plt.xlabel('wavelength (nm)')\n",
    "plt.ylabel('reflectance')\n",
    "plt.legend([str(i+1) for i in range(q)])\n",
    "plt.title('SPICE Endmembers')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# re-reval abundance maps\n",
    "P_imgs = []\n",
    "for i in range(n_em):\n",
    "    map_lin = np.zeros((n_r*n_c,))\n",
    "    map_lin[valid_array] = P[:,i]\n",
    "    P_imgs.append(np.reshape(map_lin,(n_r,n_c)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# display abundance maps\n",
    "for i in range(n_em):\n",
    "    plt.figure()\n",
    "    plt.imshow(P_imgs[i],vmin=0,vmax=1)\n",
    "    plt.colorbar()\n",
    "    plt.title('SPICE Abundance Map %d'%(i+1,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
